
rm(list=ls())
#setwd("/Users/noahbuckley/Dropbox/Cooperating/Data/EAS replication files/")
library(foreign)
library(xtable)
library(lattice)
#library(Hmisc)
library(memisc)
library(stargazer)
library(reporttools)
library(arm)
library(Hmisc)

rus <- read.dta("EAS Russia prepared data.dta")
geo <- read.dta("EAS Georgia prepared data.dta")
dta <- read.dta("EAS merged prepared data.dta")

fam <- "serif"

##########################################################################################################################
# GEORGIA ##########
##########################################################################################################################

# Rename outcome and treatment variables to something that makes more sense and matches the order that we present these experiments in the paper (they were three of eight experiments in the survey instrument, differently ordered)
geo$y_se1 <- geo$full1b
geo$y_se2 <- geo$full6
geo$y_se3 <- geo$full3
geo$y_se4 <- geo$full2b
geo$treat_se1 <- geo$q1bvar
geo$treat_se2 <- geo$q6var
geo$treat_se3 <- geo$q3var
geo$treat_se4 <- geo$q2bvar

attach(geo)

# Generate binary indicator for item non-response
q1.nr <- recode(y_se1, 1 <- NA, 0 <- c(1,2,3,4,5))
q2.nr <- recode(y_se2, 1 <- NA, 0 <- c(1,2,3,4,5))
q3.nr <- recode(y_se3, 1 <- NA, 0 <- c(1,2,3,4,5))
q4.nr <- recode(y_se4, 1 <- NA, 0 <- c(1,2,3,4,5))


#######################################################################################################################
# TABLE 1 (right-hand panels) #########################################################################################
#######################################################################################################################

tmat <- function(qdata,qvar,qnum){
	tmat.1 <- matrix(NA,6,3)
	#means and sd
	tmat.1[1,1] <- mean(qdata[qvar==1],na.rm=T)
	tmat.1[2,1] <- sd(qdata[qvar==1],na.rm=T)
	tmat.1[3,1] <- mean(qdata[qvar==2],na.rm=T)
	tmat.1[4,1] <- sd(qdata[qvar==2],na.rm=T)
	tmat.1[1,2] <- mean(qdata[qvar==3],na.rm=T)
	tmat.1[2,2] <- sd(qdata[qvar==3],na.rm=T)
	tmat.1[3,2] <- mean(qdata[qvar==4],na.rm=T)
	tmat.1[4,2] <- sd(qdata[qvar==4],na.rm=T)
	#diff in means, se
	d.13 <- t.test(qdata[qvar==1 | qvar==3] ~ qvar[qvar==1 | qvar==3])
	d.12 <- t.test(qdata[qvar==1 | qvar==2] ~ qvar[qvar==1 | qvar==2])
	d.24 <- t.test(qdata[qvar==2 | qvar==4] ~ qvar[qvar==2 | qvar==4])
	d.34 <- t.test(qdata[qvar==4 | qvar==3] ~ qvar[qvar==4 | qvar==3])
	d.14 <- t.test(qdata[qvar==4 | qvar==1] ~ qvar[qvar==4 | qvar==1])
	tmat.1[1,3] <- d.13$est[2]-d.13$est[1]
	tmat.1[2,3] <- -tmat.1[1,3]/d.13$stat
	tmat.1[5,1] <- d.12$est[2]-d.12$est[1]
	tmat.1[6,1] <- -tmat.1[5,1]/d.12$stat
	tmat.1[3,3] <- d.24$est[2]-d.24$est[1]
	tmat.1[4,3] <- -tmat.1[3,3]/d.24$stat
	tmat.1[5,2] <- d.34$est[2]-d.34$est[1]
	tmat.1[6,2] <- -tmat.1[5,2]/d.34$stat
	tmat.1[5,3] <- d.14$est[2]-d.14$est[1]
	tmat.1[6,3] <- -tmat.1[5,3]/d.14$stat
	
	tmat.1 <- matrix(tmat.1,6,3,dimnames=list(rows=c("Treatment A1","SD","Treatment A2","SD","Diff in means","SE"),cols=c("Treatment B1","Treatment B2","Diff in means")))
	tmat.out <- xtable(tmat.1,digits=2,align=c("l","c","c","c"))
	
	return(tmat.out)
}

print(tmat(y_se1, treat_se1),type="latex",file="output/geo_tout1.tex")
print(tmat(y_se2, treat_se2),type="latex",file="output/geo_tout1.tex",append=T)
print(tmat(y_se3, treat_se3),type="latex",file="output/geo_tout1.tex",append=T)
print(tmat(y_se4, treat_se4),type="latex",file="output/geo_tout1.tex",append=T)


#######################################################################################################################
# FIGURE A2 ###########################################################################################################
#######################################################################################################################

# Specify covariates we are looking for balance over; give nice labels

# age male edu material contact_any crimein12mo polstaknow often_drive working_respondent familysize urban
balvars <- cbind(age, male, edu, material, contact_any, crimein12mo, polstaknow, often_drive, working_respondent, familysize, urban)

qnames<- list(c("Age","Male","Education","Material wealth","Any contact w/ police","Crime victim in last 12mo.","Know location of police station","Drive a car often","Employed","Family size","Urban resident"),c("Question 1b","Question 2b","Question 5","Question 6"))

# Empty matrix to store all t-test p-values in
balmat.tfull.flat <- matrix(NA, dim(balvars)[2],12, dimnames=list(qnames[[1]], NULL))
for (ix in 1:dim(balvars)[2]){
	oneq <- balvars[,ix]
	balmat.tfull.flat[ix,] <- c(t.test(oneq ~ police)$p.val, t.test(oneq ~ beating)$p.val, t.test(oneq ~ policeXbeating)$p.val,  t.test(oneq ~ anonymous)$p.val, t.test(oneq ~ q6police)$p.val, t.test(oneq ~ anonXq6police)$p.val, t.test(oneq ~ duty)$p.val, t.test(oneq ~ reward)$p.val, t.test(oneq ~ dutyXreward)$p.val, t.test(oneq ~ busy)$p.val, t.test(oneq ~ highval)$p.val, t.test(oneq ~ busyXhighval)$p.val)
}

# Open pdf file, plot all p values using a variety of characters
cx <- 1.5
pdf(file="output/geobalanceplot.pdf", width=14)
par(mar=c(4,1,1,1))
plot(balmat.tfull.flat[,12], seq(1,dim(balmat.tfull.flat)[1]), pch=16, xlim=c(-0.5,1), frame.plot=F, axes=F, ylab="", xlab=expression(paste(italic("p"),"-value from ",italic("t"),"-test",sep="")), cex=cx)
points(balmat.tfull.flat[,11], seq(1,dim(balmat.tfull.flat)[1]), pch=17, cex=cx)
points(balmat.tfull.flat[,10], seq(1,dim(balmat.tfull.flat)[1]), pch=3, cex=cx)
points(balmat.tfull.flat[,9], seq(1,dim(balmat.tfull.flat)[1]), pch=8, cex=cx)
points(balmat.tfull.flat[,8], seq(1,dim(balmat.tfull.flat)[1]), pch=10, cex=cx)
points(balmat.tfull.flat[,7], seq(1,dim(balmat.tfull.flat)[1]), pch=21, cex=cx)
points(balmat.tfull.flat[,6], seq(1,dim(balmat.tfull.flat)[1]), pch=22, cex=cx)
points(balmat.tfull.flat[,5], seq(1,dim(balmat.tfull.flat)[1]), pch=23, cex=cx)
points(balmat.tfull.flat[,4], seq(1,dim(balmat.tfull.flat)[1]), pch=24, cex=cx)
points(balmat.tfull.flat[,3], seq(1,dim(balmat.tfull.flat)[1]), pch=11, cex=cx)
points(balmat.tfull.flat[,2], seq(1,dim(balmat.tfull.flat)[1]), pch=2, cex=cx)
points(balmat.tfull.flat[,1], seq(1,dim(balmat.tfull.flat)[1]), pch=13, cex=cx)
abline(v=0.05, lty=2)
abline(v=0.1, lty=3)
text(-0.5, seq(1,dim(balmat.tfull.flat)[1]), labels=rev(qnames[[1]]), cex=0.9, adj=0)
axis(1, at=seq(0,1,0.1))
legend(-0.26, 11, legend=c("Q1 Police", "Q1 Beating", "Q1 Police x Beating", "Q2 Anonymous", "Q2 Police", "Q2 Anonymous x Police", "Q3 Civic Duty", "Q3 Monetary Reward", "Q3 Civic Duty x Reward", "Q4 Six Hours Report", "Q4 High Value Theft", "Q4 Six Hours x High Value"), pch=c(16, 17, 3, 8, 10, 21, 22, 23, 24, 11, 2, 13), cex=0.9)
dev.off()



#######################################################################################################################
# FIGURES A1-A4 (right-hand plots) ####################################################################################
#######################################################################################################################

superplot.geo <- function(data,grouper,num,legendvec,cex1,horiz,vert){
	xx <- yy <- matrix(NA,4,5)
	prp.na <- num.na <- rep(NA,4)
	for (qq in 1:4) {
		h <- hist(na.omit(data[grouper==qq]), plot=FALSE)
		diffBreaks <- h$mids[2] - h$mids[1]
		tmp1 <- c(h$mids[1]-diffBreaks, h$mids, tail(h$mids,1)+diffBreaks )
		tmp2 <- c(0, h$density/2, 0)
		xx[qq,] <- round(tmp1[tmp2>0])
		yy[qq,] <- tmp2[tmp2>0]
		num.na[qq] <- length(which(is.na(data[grouper==qq])==T))
		prp.na[qq] <- num.na[qq]/(num.na[qq]+sum(h$count))
		}
		
	xx.m <- cbind(c(0,0,0,0),xx)
	means <- c(mean(na.omit(data[grouper ==1])),mean(na.omit(data[grouper ==2])),mean(na.omit(data[grouper ==3])),mean(na.omit(data[grouper ==4])))
	factor <- max(yy)/5
	means <- means*factor
	yy.m <- cbind(c(NA,NA,NA,NA),yy)
	
	par(mar=c(5,4,2,2), family=fam)  # ylim=c(-5,max(yy))
	plot(xx.m[1,], yy.m[1,], lwd=2, col="blue",type="o",pch=19,ylim=c(0,0.6),ylab="Proportion of non-missing responses",xlab="Response (\"Certainly not\" to \"Certainly\")",xaxt="n",xlim=c(-0.1,5.6),bty="n",main=
	"",lty=1,cex=1.2)
	axis(1,at=c(seq(0,5.2,1)),labels=c("DK",1,2,3,4,5))
	lines(xx.m[2,], yy.m[2,],lwd=2,col="red",type="o",pch=1,lty=2,cex=1.2)
	lines(xx.m[3,], yy.m[3,],lwd=2,col="green",type="o",pch=17,lty=3,cex=1.2)
	lines(xx.m[4,], yy.m[4,],lwd=2,col="darkgray",type="o",pch=2,lty=4,cex=1.2)
	points(rep(0,4),prp.na,pch=18,cex=1.4,col=c("blue","red","green","darkgray"))
	legend(horiz,vert,col=c("blue","red","green","darkgray"),lwd=2.3,legend=legendvec,cex=cex1,pch=c(19,1,17,2),lty=c(1,2,3,4))
	dev.copy(pdf,file=paste("output/g2013_se",num,"super1.pdf",sep=""))
	dev.off()
}

superplot.geo(y_se1, treat_se1,"1",c("G1: Policeman stealing","G2: Policeman beating","G3: police stealing","G4: police beating"),1.1,3,0.56)
superplot.geo(y_se2, treat_se2,"2",c("G1: Police + anonymity","G2: Police + no anonymity","G3: police + anonymity","G4: police + no anonymity"),1.1,2.3,0.56)
superplot.geo(y_se3, treat_se3,"4",c("G1: Low value robbery + 2hrs","G2: Low value robbery + 6hrs","G3: High value robbery + 2hrs","G4: High value robbery + 6hrs"),0.75,1.5,0.56)
superplot.geo(y_se4, treat_se4,"3",c("G1: Baseline","G2: Civic duty to report","G3: Reward of GEL5k","G4: Civic duty + reward"),1.1,2.9,0.56)



#######################################################################################################################
# TABLE X: Non-response basics and correlations ######################################################################
#######################################################################################################################

nrtable <- cor(cbind(q1.nr, q2.nr, q3.nr, q4.nr))
nrtable <- cbind(c(mean(q1.nr), mean(q2.nr), mean(q3.nr), mean(q4.nr)), nrtable)

print(xtable(nrtable,digits=3,align=c("l","c","c","c","c","c")),type="latex",file="output/geo_nonresponsebasics1.tex")


#######################################################################################################################
# TABLE X: Item response rates by demographics; table of non-response ################################################
#######################################################################################################################

# Run logit models with and without covariates
nrols.q1 <- glm(q1.nr ~ police + beating, family=binomial)
nrols.q1.2 <- glm(q1.nr ~ police + beating + age + male + edu + material + contact_any + crimein12mo + polstaknow + often_drive + working_respondent + familysize + urban, family=binomial)

nrols.q2 <- glm(q2.nr ~ duty + reward, family=binomial)
nrols.q2.2 <- glm(q2.nr ~ duty + reward + age + male + edu + material + contact_any + crimein12mo + polstaknow + often_drive + working_respondent + familysize + urban, family=binomial)

nrols.q3 <- glm(q3.nr ~ duty + reward, family=binomial)
nrols.q3.2 <- glm(q3.nr ~ duty + reward + age + male + edu + material + contact_any + crimein12mo + polstaknow + often_drive + working_respondent + familysize + urban, family=binomial)

nrols.q4 <- glm(q4.nr ~ busy + highval, family=binomial)
nrols.q4.2 <- glm(q4.nr ~ busy + highval + age + male + edu + material + contact_any + crimein12mo + polstaknow + often_drive + working_respondent + familysize + urban, family=binomial)

# Output nice latex file with these six models
labs <- c("Stranger", "Beating","Age/100","Male","Material security","Education","Interaction with police","Occurrence of crime in past 12 months","Know location of nearest police station","New to Moscow","Russian nationality", "Civic duty", "Reward", "Busy", "High value")
t1.list <- list(nrols.q1, nrols.q1.2, nrols.q2, nrols.q2.2, nrols.q3, nrols.q3.2, nrols.q4, nrols.q4.2)
stargazer(t1.list, type="latex", style="apsr", covariate.labels=NULL, font.size="scriptsize", label="nonresponsetable1", dep.var.labels=c("Question 1","Question 2","Question 3"), out="output/geo_nonresponsetable1.tex", omit.stat=c("aic", "ser" ,"f"), table.placement="h!", title="Item Non-Response by Question", column.labels="Item Non-Reponse")



#######################################################################################################################
# TABLE X: Item response rates by stratification #####################################################################
#######################################################################################################################

# Calculate mean non-response for each question in each oblast, stratum, oblastXstratum (for nice output later)
q1o <- stats:::aggregate.formula(q1.nr ~ region, data=geo, mean)

q2o <- stats:::aggregate.formula(q2.nr ~ region, data=geo, mean)

q3o <- stats:::aggregate.formula(q3.nr ~ region, data=geo, mean)

q4o <- stats:::aggregate.formula(q4.nr ~ region, data=geo, mean)

# Get sample size in each oblast and stratum
oblastsum <- stats:::aggregate.formula(rep(1, length(region)) ~ region, data=geo, sum)

# Merge all this stuff above together in the right order so we end up with one nice table
oblastfull <- merge(oblastsum, q1o, by="region")
oblastfull <- merge(oblastfull, q2o, by="region")
oblastfull <- merge(oblastfull, q3o, by="region")
oblastfull <- merge(oblastfull, q4o, by="region")
fulltable <- oblastfull

# Output this table, rename columns, print as latex using xtable
fulltable[,1] <- as.numeric(fulltable[,1])
colnames(fulltable) <- c("region", "Sampled N", "SE1 NR", "SE2 NR", "SE3 NR", "SE4 NR")
print(xtable(fulltable, digits=c(0,0,0,3,3,3,3),align=c("cllcccc")),type="latex", include.rownames=F, table.placement="htb!", size="scriptsize")


#######################################################################################################################
# TABLE X: Spillovers ################################################################################################
#######################################################################################################################

# Run six OLS models with possible combinations of outcome variables and treatment indicators from earlier experiments
# Don't need to look for spillovers with SE1 since it was first
spill.q2 <- lm(y_se2 ~ police + beating + busy + highval + duty + reward)
q2reg.spill <- lm(y_se2 ~ anonymous*q6police + police + beating + busy + highval + duty + reward)
spill.q3 <- lm(y_se3 ~ police + beating + busy + highval)
q3reg.spill <- lm(y_se3 ~ duty*reward + police + beating + busy + highval)
spill.q4 <- lm(y_se4 ~ police + beating)
q4reg.spill <- lm(y_se4 ~ busy*highval + police + beating)

# Output all six in nice latex table
labs <- c("Civic Duty (SE2)", "Reward (SE2)", "Stranger (SE1)", "Beating (SE1)", "Busy X High Value (SE1)", "Busy (SE3)", "High Value (SE3)", "Intervening Treatment 1", "Intervening Treatment 2", "Civic Duty X Reward (SE2)")
dvs <- c(rep("SE2 Outcome",4), "SE3 Outcome", "SE3 Outcome")
t2.list <- list(spill.q2, q2reg.spill, spill.q3, q3reg.spill, spill.q4, q4reg.spill)
stargazer(t2.list, type="latex", style="apsr", covariate.labels=NULL, font.size="footnotesize", label="spillovers", dep.var.labels=NULL, out="output/geo_spillovers.tex", omit.stat=c("aic", "ser" ,"f"), table.placement="h!", title="SE2 and SE3 Spillovers", notes.align="c", notes=c("Dependent variables are outcomes for Survey Experiments (SE) 2 and 3, as noted at the top of the table.", "Predictor variables are marked according to which SE they are treatments for."))



#######################################################################################################################
# TABLE X: Descriptive statistics of covariates ######################################################################
#######################################################################################################################

balvars <- cbind(age, male, edu, material, contact_any, crimein12mo, polstaknow, often_drive, working_respondent, familysize, urban)

tableContinuous(vars=balvars, stats=c("n", "min", "max", "median", "mean", "s", "na"), prec=2, col.tit=NA, cap="Covariate Descriptive Statistics", lab="descriptive_covariates", font.size="footnotesize")




detach(geo)



##########################################################################################################################
# RUSSIA ##########
##########################################################################################################################

# Rename outcome and treatment variables to something that makes more sense and matches the order that we present these experiments in the paper (they were three of eight experiments in the survey instrument, differently ordered)
rus$y_se1 <- rus$full1b
rus$y_se2 <- rus$full6
rus$y_se3 <- rus$full3
rus$y_se4 <- rus$full2b
rus$treat_se1 <- rus$q1bvar
rus$treat_se2 <- rus$q6var
rus$treat_se3 <- rus$q3var
rus$treat_se4 <- rus$q2bvar

attach(rus)

# Generate binary indicator for item non-response
q1.nr <- recode(y_se1, 1 <- NA, 0 <- c(1,2,3,4,5))
q2.nr <- recode(y_se2, 1 <- NA, 0 <- c(1,2,3,4,5))
q3.nr <- recode(y_se3, 1 <- NA, 0 <- c(1,2,3,4,5))
q4.nr <- recode(y_se4, 1 <- NA, 0 <- c(1,2,3,4,5))


#######################################################################################################################
# TABLE 1 (left-hand panels) ##########################################################################################
#######################################################################################################################

tmat <- function(qdata,qvar,qnum){
	tmat.1 <- matrix(NA,6,3)
	#means and sd
	tmat.1[1,1] <- mean(qdata[qvar==1],na.rm=T)
	tmat.1[2,1] <- sd(qdata[qvar==1],na.rm=T)
	tmat.1[3,1] <- mean(qdata[qvar==2],na.rm=T)
	tmat.1[4,1] <- sd(qdata[qvar==2],na.rm=T)
	tmat.1[1,2] <- mean(qdata[qvar==3],na.rm=T)
	tmat.1[2,2] <- sd(qdata[qvar==3],na.rm=T)
	tmat.1[3,2] <- mean(qdata[qvar==4],na.rm=T)
	tmat.1[4,2] <- sd(qdata[qvar==4],na.rm=T)
	#diff in means, se
	d.13 <- t.test(qdata[qvar==1 | qvar==3] ~ qvar[qvar==1 | qvar==3])
	d.12 <- t.test(qdata[qvar==1 | qvar==2] ~ qvar[qvar==1 | qvar==2])
	d.24 <- t.test(qdata[qvar==2 | qvar==4] ~ qvar[qvar==2 | qvar==4])
	d.34 <- t.test(qdata[qvar==4 | qvar==3] ~ qvar[qvar==4 | qvar==3])
	d.14 <- t.test(qdata[qvar==4 | qvar==1] ~ qvar[qvar==4 | qvar==1])
	tmat.1[1,3] <- d.13$est[2]-d.13$est[1]
	tmat.1[2,3] <- -tmat.1[1,3]/d.13$stat
	tmat.1[5,1] <- d.12$est[2]-d.12$est[1]
	tmat.1[6,1] <- -tmat.1[5,1]/d.12$stat
	tmat.1[3,3] <- d.24$est[2]-d.24$est[1]
	tmat.1[4,3] <- -tmat.1[3,3]/d.24$stat
	tmat.1[5,2] <- d.34$est[2]-d.34$est[1]
	tmat.1[6,2] <- -tmat.1[5,2]/d.34$stat
	tmat.1[5,3] <- d.14$est[2]-d.14$est[1]
	tmat.1[6,3] <- -tmat.1[5,3]/d.14$stat
	
	tmat.1 <- matrix(tmat.1,6,3,dimnames=list(rows=c("Treatment A1","SD","Treatment A2","SD","Diff in means","SE"),cols=c("Treatment B1","Treatment B2","Diff in means")))
	tmat.out <- xtable(tmat.1,digits=2,align=c("l","c","c","c"))
	
	return(tmat.out)
}

print(tmat(y_se1, treat_se1),type="latex",file="output/tout1.tex")
print(tmat(y_se2, treat_se2),type="latex",file="output/tout1.tex",append=T)
print(tmat(y_se3, treat_se3),type="latex",file="output/tout1.tex",append=T)
print(tmat(y_se4, treat_se4),type="latex",file="output/tout1.tex",append=T)


#######################################################################################################################
# FIGURE A1 #########################################################################################################
#######################################################################################################################

# Specify covariates we are looking for balance over; give nice labels
# age male edu material contact_any crimein12mo polstaknow often_drive working_respondent familysize urban
balvars <- cbind(age, male, edu, material, contact_any, crimein12mo, polstaknow, often_drive, working_respondent, familysize, urban)

qnames<- list(c("Age","Male","Education","Material wealth","Any contact w/ police","Crime victim in last 12mo.","Know location of police station","Drive a car often","Employed","Family size","Urban resident"),c("Question 1b","Question 2b","Question 5","Question 6"))

# Empty matrix to store all t-test p-values in
balmat.tfull.flat <- matrix(NA, dim(balvars)[2],12, dimnames=list(qnames[[1]], NULL))
for (ix in 1:dim(balvars)[2]){
	oneq <- balvars[,ix]
	balmat.tfull.flat[ix,] <- c(t.test(oneq ~ police)$p.val, t.test(oneq ~ beating)$p.val, t.test(oneq ~ policeXbeating)$p.val,  t.test(oneq ~ anonymous)$p.val, t.test(oneq ~ q6police)$p.val, t.test(oneq ~ anonXq6police)$p.val, t.test(oneq ~ duty)$p.val, t.test(oneq ~ reward)$p.val, t.test(oneq ~ dutyXreward)$p.val, t.test(oneq ~ busy)$p.val, t.test(oneq ~ highval)$p.val, t.test(oneq ~ busyXhighval)$p.val)
}

# Open pdf file, plot all p values using a variety of characters
cx <- 1.5
pdf(file="output/rusbalanceplot.pdf", width=14)
par(mar=c(4,1,1,1))
plot(balmat.tfull.flat[,12], seq(1,dim(balmat.tfull.flat)[1]), pch=16, xlim=c(-0.5,1), frame.plot=F, axes=F, ylab="", xlab=expression(paste(italic("p"),"-value from ",italic("t"),"-test",sep="")), cex=cx)
points(balmat.tfull.flat[,11], seq(1,dim(balmat.tfull.flat)[1]), pch=17, cex=cx)
points(balmat.tfull.flat[,10], seq(1,dim(balmat.tfull.flat)[1]), pch=3, cex=cx)
points(balmat.tfull.flat[,9], seq(1,dim(balmat.tfull.flat)[1]), pch=8, cex=cx)
points(balmat.tfull.flat[,8], seq(1,dim(balmat.tfull.flat)[1]), pch=10, cex=cx)
points(balmat.tfull.flat[,7], seq(1,dim(balmat.tfull.flat)[1]), pch=21, cex=cx)
points(balmat.tfull.flat[,6], seq(1,dim(balmat.tfull.flat)[1]), pch=22, cex=cx)
points(balmat.tfull.flat[,5], seq(1,dim(balmat.tfull.flat)[1]), pch=23, cex=cx)
points(balmat.tfull.flat[,4], seq(1,dim(balmat.tfull.flat)[1]), pch=24, cex=cx)
points(balmat.tfull.flat[,3], seq(1,dim(balmat.tfull.flat)[1]), pch=11, cex=cx)
points(balmat.tfull.flat[,2], seq(1,dim(balmat.tfull.flat)[1]), pch=2, cex=cx)
points(balmat.tfull.flat[,1], seq(1,dim(balmat.tfull.flat)[1]), pch=13, cex=cx)
abline(v=0.05, lty=2)
abline(v=0.1, lty=3)
text(-0.5, seq(1,dim(balmat.tfull.flat)[1]), labels=rev(qnames[[1]]), cex=0.9, adj=0)
axis(1, at=seq(0,1,0.1))
legend(-0.26, 11, legend=c("Q1 Police", "Q1 Beating", "Q1 Police x Beating", "Q2 Anonymous", "Q2 Police", "Q2 Anonymous x Police", "Q3 Civic Duty", "Q3 Monetary Reward", "Q3 Civic Duty x Reward", "Q4 Six Hours Report", "Q4 High Value Theft", "Q4 Six Hours x High Value"), pch=c(16, 17, 3, 8, 10, 21, 22, 23, 24, 11, 2, 13), cex=0.9)
dev.off()


#######################################################################################################################
# FIGURES A1-A4 (left-hand plots) #####################################################################################
#######################################################################################################################

superplot <- function(data,grouper,num,legendvec,cex1,horiz,vert){
	xx <- yy <- matrix(NA,4,5)
	prp.na <- num.na <- rep(NA,4)
	for (qq in 1:4) {
		h <- hist(na.omit(data[grouper==qq]), plot=FALSE)
		diffBreaks <- h$mids[2] - h$mids[1]
		tmp1 <- c(h$mids[1]-diffBreaks, h$mids, tail(h$mids,1)+diffBreaks )
		tmp2 <- c(0, h$density/2, 0)
		xx[qq,] <- round(tmp1[tmp2>0])
		yy[qq,] <- tmp2[tmp2>0]
		num.na[qq] <- length(which(is.na(data[grouper==qq])==T))
		prp.na[qq] <- num.na[qq]/(num.na[qq]+sum(h$count))
		}
		
	xx.m <- cbind(c(0,0,0,0),xx)
	means <- c(mean(na.omit(data[grouper ==1])),mean(na.omit(data[grouper ==2])),mean(na.omit(data[grouper ==3])),mean(na.omit(data[grouper ==4])))
	factor <- max(yy)/5
	means <- means*factor
	yy.m <- cbind(c(NA,NA,NA,NA),yy)
	
	par(mar=c(5,4,2,2), family=fam)  # ylim=c(-5,max(yy))
	plot(xx.m[1,], yy.m[1,], lwd=2, col="blue",type="o",pch=19,ylim=c(0,0.6),ylab="Proportion of non-missing responses",xlab="Response (\"Certainly not\" to \"Certainly\")",xaxt="n",xlim=c(-0.1,5.6),bty="n",main=
	"",lty=1,cex=1.2)
	axis(1,at=c(seq(0,5.2,1)),labels=c("DK",1,2,3,4,5))
	lines(xx.m[2,], yy.m[2,],lwd=2,col="red",type="o",pch=1,lty=2,cex=1.2)
	lines(xx.m[3,], yy.m[3,],lwd=2,col="green",type="o",pch=17,lty=3,cex=1.2)
	lines(xx.m[4,], yy.m[4,],lwd=2,col="darkgray",type="o",pch=2,lty=4,cex=1.2)
	points(rep(0,4),prp.na,pch=18,cex=1.4,col=c("blue","red","green","darkgray"))
	legend(horiz,vert,col=c("blue","red","green","darkgray"),lwd=2.3,legend=legendvec,cex=cex1,pch=c(19,1,17,2),lty=c(1,2,3,4))
	dev.copy(pdf,file=paste("output/r2013_se",num,"super1.pdf",sep=""))
	dev.off()
}

superplot(y_se1,treat_se1,"1",c("G1: Policeman stealing","G2: Policeman beating","G3: Stranger stealing","G4: Stranger beating"),1.1,3,0.12)
superplot(y_se2, treat_se2,"2",c("G1: Police + anonymity","G2: Police + no anonymity","G3: Stranger + anonymity","G4: Stranger + no anonymity"),1.1,2.3,0.13)
superplot(y_se3, treat_se3,"3",c("G1: Baseline","G2: Civic duty to report","G3: Reward of RUR100k","G4: Civic duty + reward"),1.1,2.9,0.13)
superplot(y_se4, treat_se4,"4",c("G1: Low value robbery + 2hrs","G2: Low value robbery + 6hrs","G3: High value robbery + 2hrs","G4: High value robbery + 6hrs"),0.75,3,0.1)


#######################################################################################################################
# TABLE X: Non-response basics and correlations ######################################################################
#######################################################################################################################

nrtable <- cor(cbind(q1.nr, q2.nr, q3.nr, q4.nr))
nrtable <- cbind(c(mean(q1.nr), mean(q2.nr), mean(q3.nr), mean(q4.nr)), nrtable)

print(xtable(nrtable,digits=3,align=c("l","c","c","c","c","c")),type="latex",file="output/rus_nonresponsebasics1.tex")


#######################################################################################################################
# TABLE X: Item response rates by demographics; table of non-response ################################################
#######################################################################################################################

# Run logit models with and without covariates
nrols.q1 <- glm(q1.nr ~ police + beating, family=binomial)
nrols.q1.2 <- glm(q1.nr ~ police + beating + age + male + edu + material + contact_any + crimein12mo + polstaknow + often_drive + working_respondent + familysize + urban, family=binomial)

nrols.q2 <- glm(q2.nr ~ duty + reward, family=binomial)
nrols.q2.2 <- glm(q2.nr ~ duty + reward + age + male + edu + material + contact_any + crimein12mo + polstaknow + often_drive + working_respondent + familysize + urban, family=binomial)

nrols.q3 <- glm(q3.nr ~ duty + reward, family=binomial)
nrols.q3.2 <- glm(q3.nr ~ duty + reward + age + male + edu + material + contact_any + crimein12mo + polstaknow + often_drive + working_respondent + familysize + urban, family=binomial)

nrols.q4 <- glm(q4.nr ~ busy + highval, family=binomial)
nrols.q4.2 <- glm(q4.nr ~ busy + highval + age + male + edu + material + contact_any + crimein12mo + polstaknow + often_drive + working_respondent + familysize + urban, family=binomial)

# Output nice latex file with these six models
labs <- c("Stranger", "Beating","Age/100","Male","Material security","Education","Interaction with police","Occurrence of crime in past 12 months","Know location of nearest police station","New to Moscow","Russian nationality", "Civic duty", "Reward", "Busy", "High value")
t3.list <- list(nrols.q1, nrols.q1.2, nrols.q2, nrols.q2.2, nrols.q3, nrols.q3.2, nrols.q4, nrols.q4.2)
stargazer(t3.list, type="latex", style="apsr", covariate.labels=NULL, font.size="scriptsize", label="nonresponsetable1", dep.var.labels=c("Question 1","Question 2","Question 3"), out="output/rus_nonresponsetable1.tex", omit.stat=c("aic", "ser" ,"f"), table.placement="h!", title="Item Non-Response by Question", column.labels="Item Non-Reponse")



#######################################################################################################################
# TABLE X: Item response rates by stratification #####################################################################
#######################################################################################################################

# Calculate mean non-response for each question in each oblast, stratum, oblastXstratum (for nice output later)
q1o <- stats:::aggregate.formula(q1.nr ~ oblast, data=rus, mean)
q1p <- stats:::aggregate.formula(q1.nr ~ stratum, data=rus, mean)
q1op <- stats:::aggregate.formula(q1.nr ~ oblast + stratum, data=rus, mean)

q2o <- stats:::aggregate.formula(q2.nr ~ oblast, data=rus, mean)
q2p <- stats:::aggregate.formula(q2.nr ~ stratum, data=rus, mean)
q2op <- stats:::aggregate.formula(q2.nr ~ oblast + stratum, data=rus, mean)

q3o <- stats:::aggregate.formula(q3.nr ~ oblast, data=rus, mean)
q3p <- stats:::aggregate.formula(q3.nr ~ stratum, data=rus, mean)
q3op <- stats:::aggregate.formula(q3.nr ~ oblast + stratum, data=rus, mean)

q4o <- stats:::aggregate.formula(q4.nr ~ oblast, data=rus, mean)
q4p <- stats:::aggregate.formula(q4.nr ~ stratum, data=rus, mean)
q4op <- stats:::aggregate.formula(q4.nr ~ oblast + stratum, data=rus, mean)

# Get sample size in each oblast and stratum
oblastsum <- stats:::aggregate.formula(rep(1, length(oblast)) ~ oblast, data=rus, sum)
stratumsum <- stats:::aggregate.formula(rep(1, length(oblast)) ~ stratum, data=rus, sum)

# Merge all this stuff above together in the right order so we end up with one nice table
stratumfull <- merge(stratumsum, q1op, by="stratum")
stratumfull <- merge(stratumfull, q2p, by="stratum")
stratumfull <- merge(stratumfull, q3p, by="stratum")
stratumfull <- merge(stratumfull, q4p, by="stratum")
oblastfull <- merge(oblastsum, q1o, by="oblast")
oblastfull <- merge(oblastfull, q2o, by="oblast")
oblastfull <- merge(oblastfull, q3o, by="oblast")
oblastfull <- merge(oblastfull, q4o, by="oblast")
fulltable <- merge(oblastfull, stratumfull, by="oblast")

# Output this table, rename columns, print as latex using xtable
fulltable[,1] <- as.numeric(fulltable[,1])
colnames(fulltable) <- c("oblast", "Sampled N", "SE1 NR", "SE2 NR", "SE3 NR", "SE4 NR", "Region ID", "Sampled N", "SE1 NR", "SE2 NR", "SE3 NR", "SE4 NR")
print(xtable(fulltable, digits=c(0,0,0,3,3,3,3,0,0,3,3,3,3),align=c("cllccccllcccc")),type="latex", include.rownames=F, table.placement="htb!", size="scriptsize")


#######################################################################################################################
# TABLE X: Spillovers ################################################################################################
#######################################################################################################################


# Question 1b	SE1		beating police/stranger
# Question 6	SE2		anonymous police
# Question 3	SE3		duty reward
# Question 2b	SE4		busy highval

# Run six OLS models with possible combinations of outcome variables and treatment indicators from earlier experiments
spill.q2 <- lm(y_se2 ~ police + beating + busy + highval + duty + reward)
q2reg.spill <- lm(y_se2 ~ anonymous*q6police + police + beating + busy + highval + duty + reward)
spill.q3 <- lm(y_se3 ~ police + beating + busy + highval)
q3reg.spill <- lm(y_se3 ~ duty*reward + police + beating + busy + highval)
spill.q4 <- lm(y_se4 ~ police + beating)
q4reg.spill <- lm(y_se4 ~ busy*highval + police + beating)

# Output all six in nice latex table
labs <- c("Civic Duty (SE2)", "Reward (SE2)", "Stranger (SE1)", "Beating (SE1)", "Busy X High Value (SE1)", "Busy (SE3)", "High Value (SE3)", "Intervening Treatment 1", "Intervening Treatment 2", "Civic Duty X Reward (SE2)")
dvs <- c(rep("SE2 Outcome",4), "SE3 Outcome", "SE3 Outcome")
t4.list <- list(spill.q2, q2reg.spill, spill.q3, q3reg.spill, spill.q4, q4reg.spill)
stargazer(t4.list, type="latex", style="apsr", covariate.labels=NULL, font.size="footnotesize", label="spillovers", dep.var.labels=NULL, out="output/rus_spillovers.tex", omit.stat=c("aic", "ser" ,"f"), table.placement="h!", title="SE2 and SE3 Spillovers", notes.align="c", notes=c("Dependent variables are outcomes for Survey Experiments (SE) 2 and 3, as noted at the top of the table.", "Predictor variables are marked according to which SE they are treatments for."))



#######################################################################################################################
# TABLE X: Descriptive statistics of covariates ######################################################################
#######################################################################################################################

balvars <- cbind(age, male, edu, material, contact_any, crimein12mo, polstaknow, often_drive, working_respondent, familysize, urban)

tableContinuous(vars=balvars, stats=c("n", "min", "max", "median", "mean", "s", "na"), prec=2, col.tit=NA, cap="Covariate Descriptive Statistics", lab="descriptive_covariates", font.size="footnotesize")


detach(rus)




##########################################################################################################################
# FIGURE 1: Merged coefficient plots #####################################################################################
##########################################################################################################################

dta$y_se1 <- dta$full1b
dta$y_se2 <- dta$full6
dta$y_se3 <- dta$full3
dta$y_se4 <- dta$full2b
dta$treat_se1 <- dta$q1bvar
dta$treat_se2 <- dta$q6var
dta$treat_se3 <- dta$q3var
dta$treat_se4 <- dta$q2bvar

attach(dta)

####

# table 1
full1 <- lm(y_se1 ~ georgia + police + beating + policeXbeating + policeXgeorgia + beatingXgeorgia + policeXbeatingXgeorgia)

# table 2
full2 <- lm(y_se2 ~ georgia + anonymous + q6police + anonXq6police + anonymousXgeorgia + q6policeXgeorgia + anonXq6policeXgeorgia)

# table 3
full3 <- lm(y_se3 ~ georgia + duty + reward + dutyXreward + dutyXgeorgia + rewardXgeorgia + dutyXrewardXgeorgia)

# table 4
full4 <- lm(y_se3 ~ georgia + busy + highval + busyXhighval + busyXgeorgia + highvalXgeorgia + busyXhighvalXgeorgia)

xsize <- 1
line.cx <- 1.5
label.cx <- 1.2

ysize <- 1.7
rus.y <- rev(seq(0.3,(xsize*ysize),(xsize*ysize/3)))
geo.y <- rev(seq(0.2,(xsize*ysize),(xsize*ysize/3)))
diff.y <- rev(seq(0.1,(xsize*ysize),(xsize*ysize/3)))

############## Figure 1a

coefs <- coef(full1)[3:8]
coef.mat <- cbind(c(coefs[1:3]),c(coefs[4:6]))
coef.dat <- data.frame(russia = coef.mat[,1], georgia = rowSums(coef.mat), diff = coef.mat[,2])
coef.dat[,3] <- -coef.dat[,3]

se <- se.coef(full1)[3:8]
se.mat <- cbind(c(se[1:3]),c(se[4:6]))
se.dat <- data.frame(russia = se.mat[,1], georgia=NA, diff = se.mat[,2])
vc <- vcov(full1)[1:8,1:8]
se.dat$georgia[1] <- sqrt(vc["policeXgeorgia","policeXgeorgia"] + vc["police","police"] + 2*vc["police","policeXgeorgia"])
se.dat$georgia[2] <- sqrt(vc["beatingXgeorgia","beatingXgeorgia"] + vc["beating","beating"] + 2*vc["beating","beatingXgeorgia"])
se.dat$georgia[3] <- sqrt(vc["policeXbeatingXgeorgia","policeXbeatingXgeorgia"] + vc["policeXbeating","policeXbeating"] + 2*vc["policeXbeating","policeXbeatingXgeorgia"])

par(mar=c(5,3,2,2), family=fam)
plot(seq(-xsize,xsize,0.1),rep(0,21),ylim=c(0,ysize),type="n",axes=0, xlab="Coefficient Estimate", ylab="", main="")
# russia
abline(v=seq(-xsize,xsize,0.1), col="lightgray", lty=3, cex=0.01)
abline(v=0, col="gray")
for (ix in 1:3){
	segments((coef.dat$russia[ix]-2*se.dat$russia[ix]),rus.y[ix],(coef.dat$russia[ix]+2*se.dat$russia[ix]),rus.y[ix],col="blue", lwd=line.cx)
	points(coef.dat$russia[ix],rus.y[ix], pch=16, col="blue", cex=1.2)
	# georgia
	segments((coef.dat$georgia[ix]-2*se.dat$georgia[ix]),geo.y[ix],(coef.dat$georgia[ix]+2*se.dat$georgia[ix]),geo.y[ix], col="red", lwd=line.cx)
	points(coef.dat$georgia[ix],geo.y[ix], pch=17, col="red", cex=1.2)
	# diff
	segments((coef.dat$diff[ix]-2*se.dat$diff[ix]),diff.y[ix],(coef.dat$diff[ix]+2*se.dat$diff[ix]),diff.y[ix], lwd=line.cx)
	points(coef.dat$diff[ix],diff.y[ix], pch=18, cex=1.4)
}
axis(1)
text(0,1.55,"Police treatment", cex=label.cx)
text(0,0.97,"Beating treatment", cex=label.cx)
text(0,0.42,"Police X Beating", cex=label.cx)
legend(-1,1,c("Russia","Georgia","R-G Difference"),col=c("blue","red","black"),lwd=2,pch=c(16,17,18), bg="white")
intcpts <- c(round(coef(full1)[1],2),round(coef(full1)[2]+coef(full1)[1],2))
intcpts.se <- c(round(se.coef(full1)[1],2),round(sqrt(vc["(Intercept)","(Intercept)"] + vc["georgia","georgia"] + 2*vc["georgia","(Intercept)"]),2))
text(0.4,1.7,paste("Control gp mean (SE): \n Russia      ",intcpts[1],"  (",intcpts.se[1],")\n"," Georgia    ",intcpts[2],"  (",intcpts.se[2],")",sep=""),cex=1,adj=c(0,1))

dev.copy(pdf,file="output/se1_coefplot.pdf")
dev.off()

############### Figure 1b

coefs <- coef(full2)[3:8]
coef.mat <- cbind(c(coefs[1:3]),c(coefs[4:6]))
coef.dat <- data.frame(russia = coef.mat[,1], georgia = rowSums(coef.mat), diff = coef.mat[,2])
coef.dat[,3] <- -coef.dat[,3]

se <- se.coef(full2)[3:8]
se.mat <- cbind(c(se[1:3]),c(se[4:6]))
se.dat <- data.frame(russia = se.mat[,1], georgia=NA, diff = se.mat[,2])
vc <- vcov(full2)[1:8,1:8]
se.dat$georgia[1] <- sqrt(vc["anonymousXgeorgia","anonymousXgeorgia"] + vc["anonymous","anonymous"] + 2*vc["anonymous","anonymousXgeorgia"])
se.dat$georgia[2] <- sqrt(vc["q6policeXgeorgia","q6policeXgeorgia"] + vc["q6police","q6police"] + 2*vc["q6police","q6policeXgeorgia"])
se.dat$georgia[3] <- sqrt(vc["anonXq6policeXgeorgia","anonXq6policeXgeorgia"] + vc["anonXq6police","anonXq6police"] + 2*vc["anonXq6police","anonXq6policeXgeorgia"])

par(mar=c(5,3,2,2), family=fam)
plot(seq(-xsize,xsize,0.1),rep(0,21),ylim=c(0,ysize),type="n", axes=0, xlab="Coefficient Estimate", ylab="", main="")
# russia
abline(v=seq(-xsize,xsize,0.1), col="lightgray", lty=3, cex=0.01)
abline(v=0, col="gray")
for (ix in 1:3){
	segments((coef.dat$russia[ix]-2*se.dat$russia[ix]),rus.y[ix],(coef.dat$russia[ix]+2*se.dat$russia[ix]),rus.y[ix],col="blue", lwd=line.cx)
	points(coef.dat$russia[ix],rus.y[ix], pch=16, col="blue", cex=1.2)
	# georgia
	segments((coef.dat$georgia[ix]-2*se.dat$georgia[ix]),geo.y[ix],(coef.dat$georgia[ix]+2*se.dat$georgia[ix]),geo.y[ix], col="red", lwd=line.cx)
	points(coef.dat$georgia[ix],geo.y[ix], pch=17, col="red", cex=1.2)
	# diff
	segments((coef.dat$diff[ix]-2*se.dat$diff[ix]),diff.y[ix],(coef.dat$diff[ix]+2*se.dat$diff[ix]),diff.y[ix], lwd=line.cx)
	points(coef.dat$diff[ix],diff.y[ix], pch=18, cex=1.4)
}# other
axis(1)
text(0,1.55,"Anonymity treatment", cex=label.cx)
text(0,0.97,"Police treatment", cex=label.cx)
text(0,0.42,"Anonymity X Police", cex=label.cx)
legend(-1,1.3,c("Russia","Georgia","R-G Difference"),col=c("blue","red","black"),lwd=2,pch=c(16,17,18), bg="white")
intcpts <- c(round(coef(full2)[1],2),round(coef(full2)[2]+coef(full2)[1],2))
intcpts.se <- c(round(se.coef(full2)[1],2),round(sqrt(vc["(Intercept)","(Intercept)"] + vc["georgia","georgia"] + 2*vc["georgia","(Intercept)"]),2))
text(0.4,1.7,paste("Control gp mean (SE): \n Russia      ",intcpts[1],"  (",intcpts.se[1],")\n"," Georgia    ",intcpts[2],"  (",intcpts.se[2],")",sep=""),cex=1,adj=c(0,1))

dev.copy(pdf,file="output/se2_coefplot.pdf")
dev.off()

############### Figure 1c

coefs <- coef(full3)[3:8]
coef.mat <- cbind(c(coefs[1:3]),c(coefs[4:6]))
coef.dat <- data.frame(russia = coef.mat[,1], georgia = rowSums(coef.mat), diff = coef.mat[,2])
coef.dat[,3] <- -coef.dat[,3]

se <- se.coef(full3)[3:8]
se.mat <- cbind(c(se[1:3]),c(se[4:6]))
se.dat <- data.frame(russia = se.mat[,1], georgia=NA, diff = se.mat[,2])
vc <- vcov(full3)[1:8,1:8]
se.dat$georgia[1] <- sqrt(vc["dutyXgeorgia","dutyXgeorgia"] + vc["duty","duty"] + 2*vc["duty","dutyXgeorgia"])
se.dat$georgia[2] <- sqrt(vc["rewardXgeorgia","rewardXgeorgia"] + vc["reward","reward"] + 2*vc["reward","rewardXgeorgia"])
se.dat$georgia[3] <- sqrt(vc["dutyXrewardXgeorgia","dutyXrewardXgeorgia"] + vc["dutyXreward","dutyXreward"] + 2*vc["dutyXreward","dutyXrewardXgeorgia"])

par(mar=c(5,3,2,2))
plot(seq(-xsize,xsize,0.1),rep(0,21),ylim=c(0,ysize),type="n",axes=0, xlab="Coefficient Estimate", ylab="", main="")
# russia
abline(v=seq(-xsize,xsize,0.1), col="lightgray", lty=3, cex=0.01)
abline(v=0, col="gray")
for (ix in 1:3){
	segments((coef.dat$russia[ix]-2*se.dat$russia[ix]),rus.y[ix],(coef.dat$russia[ix]+2*se.dat$russia[ix]),rus.y[ix],col="blue", lwd=line.cx)
	points(coef.dat$russia[ix],rus.y[ix], pch=16, col="blue", cex=1.2)
	# georgia
	segments((coef.dat$georgia[ix]-2*se.dat$georgia[ix]),geo.y[ix],(coef.dat$georgia[ix]+2*se.dat$georgia[ix]),geo.y[ix], col="red", lwd=line.cx)
	points(coef.dat$georgia[ix],geo.y[ix], pch=17, col="red", cex=1.2)
	# diff
	segments((coef.dat$diff[ix]-2*se.dat$diff[ix]),diff.y[ix],(coef.dat$diff[ix]+2*se.dat$diff[ix]),diff.y[ix], lwd=line.cx)
	points(coef.dat$diff[ix],diff.y[ix], pch=18, cex=1.4)
}# other
axis(1)
text(0,1.55,"Civic Duty Treatment", cex=label.cx)
text(0,0.97,"Reward Treatment", cex=label.cx)
text(0,0.42,"Civic Duty X Reward", cex=label.cx)
legend(-1,1.17,c("Russia","Georgia","R-G Difference"),col=c("blue","red","black"),lwd=2,pch=c(16,17,18), bg="white")
intcpts <- c(round(coef(full3)[1],2),round(coef(full3)[2]+coef(full3)[1],2))
intcpts.se <- c(round(se.coef(full3)[1],2),round(sqrt(vc["(Intercept)","(Intercept)"] + vc["georgia","georgia"] + 2*vc["georgia","(Intercept)"]),2))
text(0.4,1.7,paste("Control gp mean (SE): \n Russia      ",intcpts[1],"  (",intcpts.se[1],")\n"," Georgia    ",intcpts[2],"  (",intcpts.se[2],")",sep=""),cex=1,adj=c(0,1))

dev.copy(pdf,file="output/se3_coefplot.pdf")
dev.off()

############### Figure 1d

coefs <- coef(full4)[3:8]
coef.mat <- cbind(c(coefs[1:3]),c(coefs[4:6]))
coef.dat <- data.frame(russia = coef.mat[,1], georgia = rowSums(coef.mat), diff = coef.mat[,2])
coef.dat[,3] <- -coef.dat[,3]

se <- se.coef(full4)[3:8]
se.mat <- cbind(c(se[1:3]),c(se[4:6]))
se.dat <- data.frame(russia = se.mat[,1], georgia=NA, diff = se.mat[,2])
vc <- vcov(full4)[1:8,1:8]
se.dat$georgia[1] <- sqrt(vc["busyXgeorgia","busyXgeorgia"] + vc["busy","busy"] + 2*vc["busy","busyXgeorgia"])
se.dat$georgia[2] <- sqrt(vc["highvalXgeorgia","highvalXgeorgia"] + vc["highval","highval"] + 2*vc["highval","highvalXgeorgia"])
se.dat$georgia[3] <- sqrt(vc["busyXhighvalXgeorgia","busyXhighvalXgeorgia"] + vc["busyXhighval","busyXhighval"] + 2*vc["busyXhighval","busyXhighvalXgeorgia"])

par(mar=c(5,3,2,2), family=fam)
plot(seq(-xsize,xsize,0.1),rep(0,21),ylim=c(0,ysize),type="n",axes=0, xlab="Coefficient Estimate", ylab="", main="")
# russia
abline(v=seq(-xsize,xsize,0.1), col="lightgray", lty=3, cex=0.01)
abline(v=0, col="gray")
for (ix in 1:3){
	segments((coef.dat$russia[ix]-2*se.dat$russia[ix]),rus.y[ix],(coef.dat$russia[ix]+2*se.dat$russia[ix]),rus.y[ix],col="blue", lwd=line.cx)
	points(coef.dat$russia[ix],rus.y[ix], pch=16, col="blue", cex=1.2)
	# georgia
	segments((coef.dat$georgia[ix]-2*se.dat$georgia[ix]),geo.y[ix],(coef.dat$georgia[ix]+2*se.dat$georgia[ix]),geo.y[ix], col="red", lwd=line.cx)
	points(coef.dat$georgia[ix],geo.y[ix], pch=17, col="red", cex=1.2)
	# diff
	segments((coef.dat$diff[ix]-2*se.dat$diff[ix]),diff.y[ix],(coef.dat$diff[ix]+2*se.dat$diff[ix]),diff.y[ix], lwd=line.cx)
	points(coef.dat$diff[ix],diff.y[ix], pch=18, cex=1.4)
}# other
axis(1)
text(0,1.55,"Busy Treatment", cex=label.cx)
text(0,0.97,"High Value Theft Treatment", cex=label.cx)
text(0,0.42,"Busy X High Value", cex=label.cx)
legend(-1,1.7,c("Russia","Georgia","R-G Difference"),col=c("blue","red","black"),lwd=2,pch=c(16,17,18), bg="white")
intcpts <- c(round(coef(full4)[1],2),round(coef(full4)[2]+coef(full4)[1],2))
intcpts.se <- c(round(se.coef(full4)[1],2),round(sqrt(vc["(Intercept)","(Intercept)"] + vc["georgia","georgia"] + 2*vc["georgia","(Intercept)"]),2))
text(0.4,1.7,paste("Control gp mean (SE): \n Russia      ",intcpts[1],"    (",intcpts.se[1],")\n"," Georgia    ",intcpts[2],"  (",intcpts.se[2],")",sep=""),cex=1,adj=c(0,1))

dev.copy(pdf,file="output/se4_coefplot.pdf")
dev.off()


detach(dta)








